%% Read different versions of the same model
% by Jaromir Benes
%
% In this m-file, we create different versions of the same exogenous-growth
% model: a stationarised version (where all unit-root variables are
% transformed by dividing them by the level of productivity), and two
% unit-root (non-stationary) versions solved around different points on the
% balanced-growth path. Later on, we show that the properties of all these
% are identical. The bottom line is that we do not need to stationarised,
% or transform balanced-growth path models in any other way, to be able to
% solve them and work with them.

%% Clear workspace

clear;
close all;
home;

%% Create a baseline parameter database
%
% Choose some parameters. The parameter correspond to a model at yearly
% frequency.

P = struct();
P.g = 1.03;
P.gamma = 0.60;
P.delta = 0.10;
P.beta = 0.97;

%% Load Model Files
%
% Load both versions of the model: one version is stationarised
% (`exog_growth_stationarised.model`), the other is the same model without
% any transformation (`exog_growth_unit_root.model`), i.e. variables are in
% their original forms, not stationarised, and the model code thus
% preserves the unit root in it.

m1 = model('exog_growth_stationarised.model','assign',P) %#ok<NOPTS>
m2 = model('exog_growth_unit_root.model','assign',P) %#ok<NOPTS>

% ...
%
% Create another copy of the unit-root (non-stationary) model object for
% future use.

m3 = m2;

%% Solve for the steady states (balanced-growth paths)
%
% First, find the stationary steady state of `m1`. Then, find two different
% points on the balanced-growth path for the unit-root version of the
% model. Each of these two BGP points corresponds to a different level of
% productivity. It does not matter at all what point on the BGP is used --
% eventually, they all give exactly the same first-order solutions and the
% same results.

m1 = sstate(m1,'growth',false,'blocks',true);

m2.A = 1;
m2 = sstate(m2,'growth',true,'fixLevel',{'A'});

m3.A = 2;
m3 = sstate(m3,'growth',true,'fixLevel',{'A'});

chksstate(m1)
chksstate(m2)
chksstate(m3)

%% Check steady state (balanced-growth path)
%
% Each variable's steady state, or balanced-growth path, is described by
% two numbers: the level of the variable, and the growth rate. In BGP
% models with unit roots, there is of course no unique level toward which
% the variables would converge. What we need to know, though, is just one
% arbitrary point on the balanced-growth path. Think of it as a snapshot of
% the BGP at a particular (arbitrary) time.
%
% In IRIS, the two pieces of information that describe the steady state or
% BGP (i.e. the level and growth) are stored as complex numbers. They have,
% though, nothing to do with complex numbers -- it's just a convenient way
% of storing two things in just one number:
%
% * the real part is the level of the respective variables,
% * the imaginary part is the growth rate. For linearised variables, it is
% the difference between two consecutive periods, i.e. $\Bar x_t - \Bar
% x_{t-1}$, whereas for log-linearised variables, it is the gross rate of
% change, i.e. $\Bar x_t / \Bar x_{t-1}$.

get(m1,'sstate')
get(m2,'sstate')
get(m3,'sstate')

% ...
%
% Get only the steady-state (balanced-growth path) levels.

get(m1,'sstateLevel')
get(m2,'sstateLevel')
get(m3,'sstateLevel')

% ...
%
% Get only the steady-state (balanced-growth path) growth rates.

get(m1,'sstateGrowth')
get(m2,'sstateGrowth')
get(m3,'sstateGrowth')

%% Compute first-order dynamic solution
%
% The model `m1` is stationary, so IRIS simply computes the first-order
% Taylor expansion around the steady state, and solves for rational
% expectations. What about models `m2` and `m3`? It's surprisingly easy.
% Compute the first-order expansion around any point on the model's
% balanced-growth path. For instance, the Euler equation (with the
% expectations operator dropped for ease of notation)
%
% $$ \frac{1}{C_t} = \frac{ \beta (1+r_t) }{ C_{t+1}} $$
%
% is expanded around a point $C_t = \Bar C$, $C_{t+1} = \Bar C \cdot
% \widehat C$, and $r_t = \Bar r$, where $\Bar C$ and $\Bar r$ denote the
% respective BGP levels, and $\widehat C$ denotes the gross rate of growth.
% (these quantities are reported above). Then, solve the system for
% rational expectations.
%
% It's almost trivial to show that this procedure is equivalent to the
% following hypothetical steps:
%
% # Stationarise the model (i.e. transform the model into `m1`).
%
% # Solve the model for the transformed (stationarised) quantities, e.g.
% $y_t := Y_t/A_t$.
%
% # After solving the model, substitute back the original quantities, and
% hence obtain the law of motions for $Y_t$, etc.
%
% Instead of going through these often painful steps (recall how tedious it
% is to re-write a larger model into its stationarised equivalent), the
% original unit-root model can be handled much more easily.

m1 = solve(m1) %#ok<NOPTS>
m2 = solve(m2) %#ok<NOPTS>
m3 = solve(m3) %#ok<NOPTS>

save read_model m1 m2 m3;

%% Help on IRIS functions used in this m-file
%
% Use either `help` to display help in the command window, or `idoc`
% to display help in a HTML browser window.
%
%    help model/model
%    help model/sstate
%    help model/subsasgn
%    help model/chksstate
%    help model/get
%    help model/solve
